#|-------------------------------------------------------------------
#|Rights to this R code reserved for:                                |
#|Samuel A. Canfield                                                 |
#|Spatial Epidemiology and Ecology Research Laboratory               |
#|Department of Geography, University of Florida                     |
#|University of Florida, Gainesville, Florida, 32611 USA             |
#|  &                                                                |                                       |
#|Fredrick Tom Otieno                                                |
#|Animal Health Program, International Livestock Research Institute  |
#|P.O Box 30709 Nairobi 00100, Kenya                                 |
#---------------------------------------------------------------------
library(maps)
library(sp)
library(rgdal)
library(raster)
library(rgeos)
library(maptools)
library(mapdata)
library(gbm)
library(dismo)
library(ResourceSelection)
library(dplyr)
library(tidyverse)
library(SDMPlay)
library(spatstat)
library(xlsx)
library(MASS)
library(spatialEco)
library(pROC)
library(randomForest)
library(boot)
library(Hmisc)
library(verification)
library(pdp)
library(caret)
library(ROCR)
library(openxlsx)
library(spdep)
library(tmap)
library(GISTools)
library(aspace)
library(classInt)
library(R.utils)
library(installr)
library(devtools)
library(SDMTools)
library(GARPTools)

#Paths to SDMTools and GARPTools when needed
#install.packages("SDMTools",,"http://rforge.net/")
#install.packages("remotes")
#remotes::install_github("cghaase/GARPTools")

packages = c("maps","sp", "rgdal", "raster", "rgeos", "maptools", "mapdata", "gbm", "dismo",
             "ResourceSelection", "dplyr", "tidyverse", "SDMPlay", "spatstat", "xlsx", "MASS",
             "spatialEco", "pROC", "randomForest", "boot", "Hmisc", "verification", "pdp", "caret",
             "ROCR", "openxlsx", "spdep", "tmap", "GISTools", "aspace", "classInt", "R.utils", "installr",
             "devtools", "remotes")

## Now load or install&load all
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)

memory.limit(size=40000)

#read occurrences-presence CSV
#------------------------------
k_points_2 <- readOGR("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\lat_long_Kenya_Occurrence_201_9_3_2020.shp")
k_points_2$Species <- as.factor(k_points_2$Species)

all_kenya <- readOGR("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\arc_Full_Kenya.shp")
plot(all_kenya)

buffer_all_kenya <- readOGR("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\buffered_spatuniq_all_kenya.shp")
buffer_all_kenya <- spTransform(buffer_all_kenya,
                             crs(all_kenya))
plot(buffer_all_kenya)

#Climate Raster data Preparation
#-------------------------------------
bio4_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\bio4.tif")
bio7_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\bio7.tif")
bio13_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\bio13.tif")
llds_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\llds.tif")
pet_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\pet.tif")
pr_2_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\pr_2.tif")
pr_7_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\pr_7.tif")
pr_10_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\pr_10.tif")
pr_12_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\pr_12.tif")
slope_b <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b\\slope.tif")
plot(bio4_b)

bio4_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\bio4_rcp45_2055_ensemble.mean.v3_wc30s.tif")
bio7_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\bio7_rcp45_2055_ensemble.mean.v3_wc30s.tif")
bio13_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\bio13_rcp45_2055_ensemble.mean.v3_wc30s.tif")
llds_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\llds_rcp45_2055_ensemble.mean.v3_wc30s.tif")
pet_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\pet_rcp45_2055_ensemble.mean.v3_wc30s.tif")
pr_2_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\pr_rcp45_2055_2_ensemble.mean.v3_wc30s.tif")
pr_7_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\pr_rcp45_2055_7_ensemble.mean.v3_wc30s.tif")
pr_10_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\pr_rcp45_2055_10_ensemble.mean.v3_wc30s.tif")
pr_12_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\pr_rcp45_2055_12_ensemble.mean.v3_wc30s.tif")
slope_e45 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45\\slope.tif")
plot(pr_10_e45)

bio4_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\bio4_rcp85_2055_ensemble.mean.v3_wc30s.tif")
bio7_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\bio7_rcp85_2055_ensemble.mean.v3_wc30s.tif")
bio13_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\bio13_rcp85_2055_ensemble.mean.v3_wc30s.tif")
llds_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\llds_rcp85_2055_ensemble.mean.v3_wc30s.tif")
pet_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\pet_rcp85_2055_ensemble.mean.v3_wc30s.tif")
pr_2_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\pr_rcp85_2055_2_ensemble.mean.v3_wc30s.tif")
pr_7_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\pr_rcp85_2055_7_ensemble.mean.v3_wc30s.tif")
pr_10_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\pr_rcp85_2055_10_ensemble.mean.v3_wc30s.tif")
pr_12_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\pr_rcp85_2055_12_ensemble.mean.v3_wc30s.tif")
slope_e85 <- raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85\\slope.tif")
plot(pr_10_e85)

ma_bio4_e45 <- raster::mask(bio4_e45, all_kenya)
plot(ma_bio4_e45)
cr_bio4_e45 <- crop(ma_bio4_e45, extent(all_kenya))
plot(cr_bio4_e45)

cr_bio4_b <-raster:: resample(bio4_b, cr_bio4_e45)
cr_bio7_b <- raster::resample(bio7_b, cr_bio4_e45)
cr_bio13_b <- raster::resample(bio13_b, cr_bio4_e45)
cr_llds_b <- raster::resample(llds_b, cr_bio4_e45)
cr_pet_b <- raster::resample(pet_b, cr_bio4_e45)
cr_pr_2_b <- raster::resample(pr_2_b, cr_bio4_e45)
cr_pr_7_b <- raster::resample(pr_7_b, cr_bio4_e45)
cr_pr_10_b <- raster::resample(pr_10_b, cr_bio4_e45)
cr_pr_12_b <- raster::resample(pr_12_b, cr_bio4_e45)
cr_slope_b <- raster::resample(slope_b, cr_bio4_e45)

ma_bio7_e45 <- raster::mask(bio7_e45, all_kenya)
cr_bio7_e45 <- raster::resample(ma_bio7_e45, cr_bio4_e45)
ma_bio13_e45 <- raster::mask(bio13_e45, all_kenya)
cr_bio13_e45 <- raster::resample(ma_bio13_e45, cr_bio4_e45)
ma_llds_e45 <- raster::mask(llds_e45, all_kenya)
cr_llds_e45 <- raster::resample(ma_llds_e45, cr_bio4_e45)
ma_pet_e45 <- raster::mask(pet_e45, all_kenya)
cr_pet_e45 <- raster::resample(ma_pet_e45, cr_bio4_e45)
ma_pr_2_e45 <- raster::mask(pr_2_e45, all_kenya)
cr_pr_2_e45 <- raster::resample(ma_pr_2_e45, cr_bio4_e45)
ma_pr_7_e45 <- raster::mask(pr_7_e45, all_kenya)
cr_pr_7_e45 <- raster::resample(ma_pr_7_e45, cr_bio4_e45)
ma_pr_10_e45 <- raster::mask(pr_10_e45, all_kenya)
cr_pr_10_e45 <- raster::resample(ma_pr_10_e45, cr_bio4_e45)
ma_pr_12_e45 <- raster::mask(pr_12_e45, all_kenya)
cr_pr_12_e45 <- raster::resample(ma_pr_12_e45, cr_bio4_e45)
ma_slope_e45 <- raster::mask(slope_e45, all_kenya)
cr_slope_e45 <- raster::resample(ma_slope_e45, cr_bio4_e45)

ma_bio4_e85 <- raster::mask(bio4_e85, all_kenya)
cr_bio4_e85 <- raster::resample(ma_bio4_e85, cr_bio4_e45)
ma_bio7_e85 <- raster::mask(bio7_e85, all_kenya)
cr_bio7_e85 <- raster::resample(ma_bio7_e85, cr_bio4_e45)
ma_bio13_e85 <- raster::mask(bio13_e85, all_kenya)
cr_bio13_e85 <- raster::resample(ma_bio13_e85, cr_bio4_e45)
ma_llds_e85 <- raster::mask(llds_e85, all_kenya)
cr_llds_e85 <- raster::resample(ma_llds_e85, cr_bio4_e45)
ma_pet_e85 <- raster::mask(pet_e85, all_kenya)
cr_pet_e85 <- raster::resample(ma_pet_e85, cr_bio4_e45)
ma_pr_2_e85 <- raster::mask(pr_2_e85, all_kenya)
cr_pr_2_e85 <- raster::resample(ma_pr_2_e85, cr_bio4_e45)
ma_pr_7_e85 <- raster::mask(pr_7_e85, all_kenya)
cr_pr_7_e85 <- raster::resample(ma_pr_7_e85, cr_bio4_e45)
ma_pr_10_e85 <- raster::mask(pr_10_e85, all_kenya)
cr_pr_10_e85 <- raster::resample(ma_pr_10_e85, cr_bio4_e45)
ma_pr_12_e85 <- raster::mask(pr_12_e85, all_kenya)
cr_pr_12_e85 <- raster::resample(ma_pr_12_e85, cr_bio4_e45)
ma_slope_e85 <- raster::mask(slope_e85, all_kenya)
cr_slope_e85 <- raster::resample(ma_slope_e85, cr_bio4_e45)

predictors_anthrax_b <-stack(cr_bio4_b, cr_bio7_b, cr_bio13_b, cr_llds_b, cr_pet_b, cr_pr_2_b, cr_pr_7_b, cr_pr_10_b, cr_pr_12_b, cr_slope_b)
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b")
writeRaster(predictors_anthrax_b, filename=names(predictors_anthrax_b), bylayer=TRUE,format="GTiff",overwrite=TRUE)

predictors_anthrax_e45 <- stack(cr_bio4_e45, cr_bio7_e45, cr_bio13_e45, cr_llds_e45, cr_pet_e45, cr_pr_2_e45, cr_pr_7_e45, cr_pr_10_e45, cr_pr_12_e45,
                                cr_slope_e45)
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45")
writeRaster(predictors_anthrax_e45, filename=names(predictors_anthrax_e45), bylayer=TRUE,format="GTiff",overwrite=TRUE)

predictors_anthrax_e85 <- stack(cr_bio4_e85, cr_bio7_e85, cr_bio13_e85, cr_llds_e85, cr_pet_e85, cr_pr_2_e85, cr_pr_7_e85, cr_pr_10_e85, cr_pr_12_e85,
                                cr_slope_e85)
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85")
writeRaster(predictors_anthrax_e85, filename=names(predictors_anthrax_e85), bylayer=TRUE,format="GTiff",overwrite=TRUE)


#Call climate covariate Rasters------------------------------

setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_b") #Tiffs path
tifFiles <- Sys.glob ('*.tif')
predictors_anthrax_b <-stack(tifFiles,quick
                           =F)
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e45") #Tiffs path
tifFiles <- Sys.glob ('*.tif')
predictors_anthrax_e45 <-stack(tifFiles,quick                         =F)

setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\Covariates_e85") #Tiffs path
tifFiles <- Sys.glob ('*.tif')
predictors_anthrax_e85 <-stack(tifFiles,quick
                           =F)
                          
#Prepare Spatially Uniqupoint data for model runs: All of Kenya
#Determine centroid locations of raster cells with presence locations to remove duplicate presence points
#----------------------------------------------------------------------------------------------------------
su_bio13_b <- predictors_anthrax_b@layers[[1]]
spatuniq <-  centroid(x = su_bio13_b, points = k_points_2, species = k_points_2$Species)
plot(spatuniq)
spatuniq_2 <- spatuniq
spatuniq_2 <- data.frame(spatuniq@data[["Longitude"]], spatuniq@data[["Latitude"]], spatuniq@data[["Species"]])
colnames(spatuniq_2)<- c("X","Y", "Species")
spatuniq_3 <- SpatialPointsDataFrame(spatuniq_2[,c("X", "Y")], spatuniq_2[,1:3])
writeOGR(spatuniq_3, layer = "spatuniq_Kenya_Occurrences", "C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\spatuniq_Kenya_Occurrences.shp", driver="ESRI Shapefile")

#Kenya Baseline------------------------

#loop run steps
set.seed(100)
mylist_b<-list()
ivlist_b <- list()
testlist_b <- list()
trainlist_b <- list()
qs.list_b <- list()
nrep= 2
for (r in 1:nrep){
  setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base")
  nam <- print(paste("run", r))
  
  #generate random points
  #-------------------------
  anthrax_pseudos<-spsample(buffer_all_kenya, n=201,"random")
  anthrax_pseudos_df<-data.frame(anthrax_pseudos)
  anthrax_pseudos_df<-anthrax_pseudos_df[,1:2]
  points(anthrax_pseudos)
  
  #training(75%) and Testing (25%) 
  #---------------------------------
  #Presence data
  trainIndex <- createDataPartition(spatuniq_2$Species, p=0.75, list = FALSE)
  presence_Train <- spatuniq_2[trainIndex,]
  presence_Test <- spatuniq_2[-trainIndex,]
  presence_Train$Species <- NULL
  presence_Test$Species <- NULL
  
  #absence data
  anthrax_pseudos_df$species <- c("anthrax")
  trainIndex_a <- createDataPartition(anthrax_pseudos_df$species, p=0.75, list = FALSE)
  abs_Train <- anthrax_pseudos_df[trainIndex_a,]
  abs_Test <- anthrax_pseudos_df[-trainIndex_a,]
  abs_Train$species <- NULL
  abs_Test$species <- NULL
  
  #extracting train presence and absence data
  #-------------------------------------------
  file1<-raster::extract(predictors_anthrax_b,presence_Train,method='simple')
  outcome<-rep(1,dim(file1)[[1]])
  file1<-cbind(outcome,file1)
  
  file2<-raster::extract(predictors_anthrax_b,abs_Train,method='simple')
  outcome<-rep(0,dim(file2)[[1]])
  file2<-cbind(outcome,file2)
  
  train_data<-as.data.frame(rbind(file1,file2))
  train_data<-na.omit(train_data)
  train_data2<-train_data
  colnames(train_data2) <- c('outcome','bio13','bio4','bio7','llds','pet','pr_10','pr_12','pr_2','pr_7','slope')
  getwd()
  write.csv(train_data2,paste0("train","_", r, ".csv"),row.names=TRUE,col.names=FALSE)
  trainlist_b[[r]] <- train_data2
  
  #extracting test presence and absence data
  #----------------------------------------------
  pre_testdata<-raster::extract(predictors_anthrax_b,presence_Test,method='simple')
  outcome<-rep(1,dim(pre_testdata)[[1]])
  pre_testdata<-cbind(outcome,pre_testdata)
  
  abs_testdata<-raster::extract(predictors_anthrax_b,abs_Test,method='simple')
  outcome<-rep(0,dim(abs_testdata)[[1]])
  abs_testdata<-cbind(outcome,abs_testdata)
  
  test_data<-as.data.frame(rbind(pre_testdata,abs_testdata))
  test_data<-na.omit(test_data)
  test_data2 <- test_data
  colnames(test_data2) <- c('outcome','bio13','bio4','bio7','llds','pet','pr_10','pr_12','pr_2','pr_7','slope')
  write.csv(test_data2,paste0("test","_", r, ".csv"),row.names=TRUE,col.names=FALSE)
  testlist_b[[r]] <- test_data2
  
  #fitting model and selecting variables
  #######################################
  col<-ncol(train_data)
  
  brt_step2 <-gbm.step(data = train_data2, gbm.x = c(2:col), gbm.y = 1,family = "bernoulli", tree.complexity = 5,learning.rate = 0.001, bag.fraction = 0.5,max.trees = 2500,n.folds = 10)
  
  #obtaining and plotting variable influence
  #-----------------------------------------
  ivlist_b[[r]]<-data.frame(summary(brt_step2))
  
  #prediction 
  ###################

  #storing brt objects
  #--------------------
  brt_step.bs <- try(gbm.step(data=train_data2[sample(NROW(train_data), NROW(train_data), replace=T),], gbm.x =c(2:col),
                              gbm.y = 1, family = "bernoulli", tree.complexity = 5, learning.rate = 0.001, bag.fraction = 0.5,max.trees = 2500,n.folds = 10,
                              verbose=TRUE, silent=FALSE, plot.main=TRUE))
  if (class(brt_step.bs) == "try-error") next
  qs.list_b[[r]] <- brt_step.bs
  cat("This is replicate number ", r, "\n")
  save(qs.list_b, file="qs.list.Rdata")
  rm(brt_step.bs)
  
  #storing run AUCs
  #----------------
  AUC<-brt_step2$cv.statistics$discrimination.mean
  mylist_b<-append(mylist_b,AUC)
  save(mylist_b, file="mylist.Rdata")
  save(ivlist_b, file="ivlist.Rdata")
  save(testlist_b, file="testlist.Rdata")
  save(trainlist_b, file="trainlist.Rdata")
}
#create CV AUC excel file
#--------------------------

write.xlsx(unlist(mylist_b),"C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\AUC.xlsx")

#Test_Data_Evaluation-------------------
#Model evaluation using testdata 
evallist_b<- list()
auclist_b <- list()
brt_p_2list_b <- list()
brt_predlist_b <- list()
for (r in 1:length(qs.list_b)) {
  print(paste("testpred", r))
  
  brt_prediction <- predict(qs.list_b[[r]], testlist_b[[r]], n.trees = qs.list_b[[r]]["n.trees"], type = "response")
  calc.deviance(obs=testlist_b[[r]]$outcome, pred=brt_prediction, calc.mean=TRUE)
  eval_data <- cbind(testlist_b[[r]]$outcome, brt_prediction)
  presences <- eval_data[eval_data[,1]==1, 2]
  abscenses <- eval_data[eval_data[,1]==0, 2]
  evaluation <- dismo::evaluate(p=presences, a=abscenses)
  evaluation
  t_AUC <- evaluation@auc
  auclist_b <- append(auclist_b, t_AUC)
  evallist_b[[r]] <- evaluation
  brt_predlist_b[[r]] <- brt_prediction
  brt_p_2 <- data.frame(brt_prediction)
  brt_p_2list_b[[r]] <- brt_p_2
  save(evallist_b, file="C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\evallist.Rdata")
  save(auclist_b, file="C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\auclist.Rdata")
  save(brt_p_2list_b, file="C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\brt_p_2list.Rdata")
  save(brt_predlist_b, file="C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\brt_predlist.Rdata")
}

write.xlsx(unlist(auclist_b), ("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\test_AUC.xlsx"))
#Kenya Baseline Landscape Prediction-----------------------------------------

#Training Runs to predict onto the landscape
b_ptlist <- list()
for (r in 1:length(qs.list_b)) {
  setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base\\Landscape_Predictions_B\\Rasters")
  print(paste("pred", r))
  
  b_pt <- predict(predictors_anthrax_b, qs.list_b[[r]], n.trees = qs.list_b[[r]]["n.trees"], type = "response")
  myfile <- paste0("pred","_", r, ".tif")
  writeRaster(b_pt, filename = myfile, format = "GTiff", overwrite = TRUE)
  b_ptlist[[r]] <- b_pt
}

setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_base")
save(b_ptlist, file="b_ptlist.Rdata")
#Mean of Predictions
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Rasters")
tifFiles <- Sys.glob('*.tif') 
predicted<-stack(tifFiles, quick=F)
mean_narm = function(x,...){mean(x,na.rm=TRUE)}
Predict_mean<- do.call(overlay, c(predicted, fun = mean_narm))
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI")
writeRaster(Predict_mean,filename='mean_e85', format="GTiff", overwrite=TRUE)
meanPred_e85<-raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI\\mean_e85.tif")
plot(meanPred_e85)

#Upper 97.5% of predictions
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Rasters")
tifFiles <- Sys.glob('*.tif') 
predicted<-stack(tifFiles, quick=F)
uci_narm = function(x,...){quantile(x, probs = c(0.975),na.rm=TRUE)}
Predict_uci<- do.call(overlay, c(predicted, fun = uci_narm))
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI")
writeRaster(Predict_uci,filename='uci_e85', format="GTiff", overwrite=TRUE)
uciPred_e85<-raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI\\uci_e85.tif")
plot(uciPred_e85)

#Lower 2.5% of predictions
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Rasters")
tifFiles <- Sys.glob('*.tif') 
predicted<-stack(tifFiles, quick=F)
lci_narm = function(x,...){quantile(x, probs = c(0.025),na.rm=TRUE)}
Predict_lci<- do.call(overlay, c(predicted, fun = lci_narm))
setwd("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI")
writeRaster(Predict_lci,filename='lci_e85', format="GTiff", overwrite=TRUE)
lciPred_e85<-raster("C:\\RF_data\\Inputs\\Kenya_data\\ClimateScenarios\\Climate_Kenya_sc_v3\\train_85\\Landscape_Predictions_85\\Averaged_UCI_LCI\\lci_e85.tif")
plot(lciPred_e85)







